##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(parallel)
library(dplyr)

##########################################################################################
option_list <- list(
    make_option(c("--rna_file"), type = "character"),
    make_option(c("--gene_list_file"), type = "character"),
    make_option(c("--cluster"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 单细胞表达文件
    rna_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined.annotationCellType.qc.Rdata"

    ## 感兴趣的基因列表
    gene_list_file <- "~/20231121_singleMuti/results/tf_regulators/gene_list/cluster5.gene_interest.list"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/expression_correlation"

    ## 细胞类型
    cluster <- "cluster5"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_file <- opt$rna_file
cluster <- opt$cluster
out_path <- opt$out_path
gene_list_file <- opt$gene_list_file

dir.create( out_path , recursive = T)

###########################################################################################

a <- load(rna_file)
DefaultAssay(scrnat) <- "RNA"
## scrnat

## 需要计算的基因列表
gene_list <- fread(gene_list_file , header = F)$V1

###########################################################################################
## 提取感兴趣的细胞
cluster_use <- gsub( "cluster" , "" , cluster )
cell_use <- names(scrnat$cell[which(scrnat$seurat_clusters==cluster_use)])

## 提取表达矩阵
exprSet <- GetAssayData(object=scrnat,assay="RNA",layer="data")
exprSet <- as.matrix(exprSet)
## 特异细胞类型和特定基因集合
exprSet <- exprSet[gene_list,cell_use]

cell_type <- unique(scrnat$cell_type[cell_use])

## 移除变量
rm(scrnat)

###########################################################################################
## 计算基因间表达相关性
cluster_cor <- bind_rows(mclapply(1:(nrow(exprSet)-1),function(i){
    print(i)
    tmp_res <- c()
    for( j in (i+1):nrow(exprSet)){

        geneA <- rownames(exprSet)[i]
        geneB <- rownames(exprSet)[j]

        cor_res <- cor.test(exprSet[i,],exprSet[j,],method = 'spearman')
        tmp_res <- bind_rows( tmp_res , data.frame(cell_type = cell_type , cluster = cluster , geneA = geneA , geneB = geneB , correlation=cor_res$estimate,pvalue=cor_res$p.value) )
    }
    return(tmp_res)
},mc.cores=5))

out_file <- paste0(out_path , "/" , cluster , "_expressionCorrelation.tsv")
write.table(cluster_cor , out_file , row.names = F , sep = "\t" , quote = F )